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Abstract. We investigate the evolution of rigidly and differentially rotating protoneutron stars during the first twenty seconds 
of their life. We solve the equations describing stationary axisymmetric configurations in general relativity coupled to a finite 
temperature, relativistic equation of state, to obtain a sequence of quasi-equilibrium configurations describing the evolution of 
newly born neutron stars. The initial rotation profiles have been taken to mimic the situation found immediately following the 
gravitational collapse of rotating stellar cores. By analyzing the output of several models, we estimate that the scale of variation 
of the angular velocity in a newly born neutron star is of the order of 7-10 km. We obtain the maximum rotation frequency that 
can be reached as the protoneutron stars deleptonizes and cools down, as well as other relevant parameters such as total angular 
momentum or the instability parameter \T/W\. Our study shows that imposing physical constraints (conservation of baryonic 
mass and angular momentum) and choosing reasonable thermodynamic al profiles as the star evolves gives results consistent 
with the energetics of more complex simulations of non-rotating protoneutron stars. It appears to be unlikely that newly born 
protoneutron stars formed in nearly axisymmetric core collapses reach the critical angular velocity to undergo the bar mode 
instability. They could, however, undergo secular or low \T/W\ rotational instabilities a few seconds after birth, resulting in a 
strong emission of gravitational waves retarded with respect to the neutrino luminosity peak. We also found that the geometry 
of strongly differentially rotating protoneutron stars can become toroidal-like for large values of the angular velocity, before 
reaching the mass shedding limit. 
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1. Introduction 

The first minute of life of non-rotating protoneutron stars 
(PNS) has been studied in detail since the mid 80s by several 
authors (Burrows & Lattimer 1986; Keil & Janka 1986; Pons 
et al. 1999, 2000, 2001). The initial hot, lepton rich remnant, 
left behind following a successful core collapse supernova is 
known to evolve in a timescale of tens of seconds to form a 
cold (T < 10 10 K) catalyzed neutron star. This process is usu- 
ally classified in three main stages: i) the mantle contraction 
during which fast cooling of the outer regions takes place in 
about 0.5 s, with probably significant accretion; ii) deleptoniza- 
tion and consequently heating of the internal core as energetic 
neutrinos diffuse out leaving most of their chemical energy 
(analogously to Joule effect) on the way out; and iii) cooling 
by means of diffusion of (mostly) thermal neutrinos, resulting 
in a decrease of temperature from about 40-50 MeV to below 
2-4 MeV, point at which the star becomes transparent to neutri- 
nos. 
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Less is known about the early evolution of PNSs that 
are rapidly rotating, because a fully consistent study is a 
formidable task. It requires solving the neutrino transport equa- 
tions in at least two-dimensions, coupled to general relativistic 
hydrodynamics that describes the fluid motion in rotating rel- 
ativistic stars. In addition, the microphysical inputs also need 
to be treated carefully. As has been shown previously (Pons et 
al. 1999), internal consistency between the equation of state of 
dense matter and the neutrino opacities describing the interac- 
tion between neutrinos and matter is needed to achieve reliable 
results. 

A few attempts to treat the problem in a simplified way ex- 
ist (Romero et al. 1992; Goussard et al. 1997, 1998; Sumiyoshi 
et al. 1999; Strobel et al. 1999; Yuan & Heyl 2003) based 
on the analysis of temporal sequences of quasi-equilibrium 
models of PNS, imposing some ad-hoc thermodynamical pro- 
files {i.e., constant temperature, constant entropy and/or con- 
stant neutrino fraction). The assumption of quasi-equilibrium 
is well justified, since the hydrodynamical timescale (10~ 3 s) is 
much smaller than the timescale in which substantial thermo- 
dynamical changes occur (diffusion timescale « 1 s). However, 
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it is unclear whether or not an isothermal or isentropic pro- 
file can effectively mimic reality. In addition, there are some 
physical constraints such as conservation of energy (the total 
gravitational mass must decrease consistently with the neu- 
trino luminosity), or conservation of angular momentum, that 
need to be satisfied when one studies a temporal sequence of 
quasi-equilibrium models of PNS. These constraints help to re- 
duce the parameter space as one calculates the thermodynam- 
ical evolution, although there is much uncertainty about the 
transport of angular momentum. In general angular momentum 
losses by neutrinos are supposed to be small (Kazana 1977). 
Neutrinos could redistribute or take away a fraction of the ini- 
tial angular momentum, but quantifying this amount needs of 
multidimensional transport simulations. It is equally controver- 
sial how important is turbulent transport, or the role of gravita- 
tional radiation. Because of these many unresolved issues it is 
usually assumed that the angular momentum is approximately 
conserved during the Kelvin-Helmholtz stage. In addition, the 
presence of large magnetic field could be crucial. It has been re- 
cently shown that turbulent mean-field dynamo action can be 
effective for PNSs with periods shorter than 1 s, which would 
generate very strong magnetic fields in the interior (Bonanno et 
al. 2003). 

The whole problem is far to be solved, and will require 
some serious computational effort but, in the meantime, the 
intention of this paper is to begin to understand qualitatively 
what kind of outcome will be revealed by such numerical sim- 
ulations, prior to engaging large scale simulations. In this line, 
we improve in several respects the few existing previous works. 
First, instead of isentropic or isothermal models, we use realis- 
tic profiles coming from ID simulations, conveniently rescal- 
ing the temperature and chemical profiles as functions of den- 
sity; second, we do not restrict ourselves to the rigid rotation 
case, as most of previous studies, and include differential ro- 
tation in our analysis; last but not least, we check that fixing 
the baryonic mass and the angular momentum results in some 
reasonable time dependence to the integrated relevant physical 
quantities such as the gravitational mass. In this way, we can 
calculate a complete evolutionary sequence of a PNS . 

The plan of the paper is as follows. In 2 we briefly review 
the relativistic description of rotating stars, and the numerical 
method used to obtain equilibrium configurations. A similar 
brief description of the thermodynamics and the equation of 
state is presented in 3, since for both issues the details have 
been published elsewhere. Section 4 discusses the choice of 
parameters for the initial models. Our results are presented in 
section 5. In the first part we describe rigidly rotating PNSs, 
while the description of differentially rotating PNSs is given in 
5.2. In both cases, we give all the relevant parameters for max- 
imally rotating configurations as well as the evolution of more 
realistic sequences. Finally, our main conclusions and findings 
are summarized and discussed in 6. 

2. Numerical code for stationary axisymmetric 
spacetime 

Since the time scale for the matter and the gravitational field to 
settle in stationary equilibrium is much shorter than all other 



involved times 1 , a simplifying but plausible approximation 
consists of mimicking the time evolution of the spacetime 
by obtaining a collection of stationary, axisymmetric and 
asymptotically flat snapshots, all of them being solutions of 
the Einstein equations, but with a source term that depends on 
the snapshot. 

This approximation has already been adopted by several 
authors in the past (Goussard et al. 1997, 1998; Strobel et 
al. 1999; Yuan & Heyl 2003), our code being moreover very 
similar to the one of Goussard et al. (1997, 1998), since we 
also employ a fully relativistic spectral code based on the work 
of Bonazzola et al. (1993). In the following, this article will be 
referred to as BGSM, and we address the interested reader to 
it for more details about the algorithms. Let us just remind that 
the BGSM approach uses the 3 + 1 formulation of general rel- 
ativity, maximal slicing-quasi-isotropic coordinates, with the 
assumption that matter is free of convective motions. The latter 
hypothesis is probably too strong, since PNSs are known to be 
convectively unstable (Keil et al. 1996; Miralles et al. 2002), 
but it implies the circularity of the stationary axisymmetric 
spacetime (Carter 1969) and greatly simplifies the calculations. 

Under these approximations, the metric is written 2 

ds 2 = g tlv dx ti Ax v 

= -[N 2 - N^N^dt 2 -IN^dtdtp 

A 4 

+ (dr 2 + r 2 d® 2 + B 4 r 2 sin 2 $ dip 2 ) , (1) 

where we have introduced the usual notation of the 3 + 1 
formalism: N is the lapse, is the third component of the 
shift 3-vector, with N v =g vi N l = A 4 B 2 N v r 2 sin 2 # being the 
covariant ^-component of the latter. Thanks to the hypothesis 
of axisymmetry and stationarity, the functions N,N' P ,A and B 
depend only on r and § . 

Despite the general formalism and structure of the code is 
very similar to that of Goussard et al. (1997, 1998), there are a 
few differences that should be mentioned: 

- we use a more recent version of the spectral BGSM code 
(see Gourgoulhon et al. 1999), in which at least two shells 
map the interior of the star, in order to better describe the 
non-spherical surface. This code is based on the C++ li- 
brary LOPvENE 3 , a software package for numerical relativ- 
ity freely available under GNU license; 

- the equation of state (EOS) is written in tabular form and 
interpolated using cubic Hermite polynomials, as defined 
in Nozawa et al. (1998), according to a general technique 
presented in Swesty (1996) that preserves thermodynami- 
cal consistency (see also further remarks below); 

- during the evolution of PNSs we do not restrict ourselves to 
3 or 4 typical evolutionary snapshots that are just isentropic 

1 unless a gravitational dynamical instability occurs. 

2 For all this article, we are using natural units in which the speed 
of light is equal to unity. 

3 http://www.lorene.obspm.fr 
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or isothermal profiles. Instead we adopt more realistic pro- 
files, rescaled from the results of ID simulations with neu- 
trino transport (Pons et al. 1999, more details are given in 
the next section). 

For clarity, let us be more specific about the last remark. 
From the hydrodynamical point of view, one must deal with 
the conditions for the existence of stationary solutions of the 
Einstein equations in presence of a rotating perfect fluid (see 
also BGSM and Goussard et al.). The main issue is that the rel- 
ativistic Euler equations, derived from the conservation of the 
energy-momentum tensor projected into a space-like 3-surface, 
are reduced in the case of stationary rotating motion to 



djP 
e + P 

where 



+ di In 



FdiCl, 



(2) 



- P is the pressure and e the total energy density, both mea- 
sured in the comoving frame; 

- N is the already defined lapse function, T the Lorentz factor 
linking the Eulerian observer and the comoving observer, 
n^u^/u' the angular velocity, and F = 11^11' with u being 
the 4-velocity of the fluid; 

- the index i is either 1 or 2, since all the previous quantities 
depend only on r and §. 

As explained in BGSM Eq. (2) does not admit a first inte- 
gral for arbitrary thermodynamics or rotation profiles Q[r, ■&], 
but it does in some simple cases. 

Concerning rotation, the right hand side of equation (2) is a 
first integral if either Q is constant in space (rigid rotation), or 
the function F can be locally written as F = F[Q.]. The latter 
leads to the relation 



F[£2] 



A 4 B 2 r 2 sin 2 (Q - N*) 
N 2 - A 4 B 2 r 2 sin 2 (Q - N*) 2 



0, 



(3) 



between O and the coefficients of the metric, which allows the 
determination of the profile of angular velocity. 

For the differential rotation, we adopt the law proposed by 
Komatsu et al. (1989) and used by many authors (Komatsu et 
al. 1989; Goussard et al. 1997, 1998; Baumgarte et al. 2000) 



F[£2] 



Q), 



(4) 



where Q. c is the limit of the function Q[r, -&] on the rotation axis 
and is a parameter with the dimension of the square of a 
length, noted A in Komatsu et al. (1989). We shall discuss later 
our choices for the values of these parameters. It must be re- 
marked that despite its simplicity, this law gives reasonable ro- 
tation profiles qualitatively similar to those obtained from core 
collapse simulations (see section 4). 

Regarding the EOS, it can be shown that isentropic or 
isothermal profiles are sufficient conditions for Eq. (2) to be 
integrable. Yet another possibility which is less restrictive ex- 
ists : assuming that the total energy density e can be written as 
an effective function of the pressure e = e[P]. Indeed, in this 
case, even if the actual EOS is not barotropic, assuming that the 



EOS is a one-parameter function ensures the existence of the 
first integral 



H + In 



F[Q] dQ = const. 



where the function H is defined by 



FL-- 



! 



dP 



e[P] + P ' 



(5) 



(6) 



with suitable boundary conditions. In this work, we have 
rescaled the temperature and composition profiles assuming, at 
each evolutionary time, that they depend only on the baryon 
density n B . In this way, the EOS behaves effectively as a 
barotrope, and the pseudo-enthalpy function (6) can be written 
as 



r 

H = 

Jo 



1 



dP[n B ] 

e[n B ] + P[n B ] 



dn B , 



(7) 



with the effective functions P[n B ] and e[n B ] depending only on 
the baryonic density n B . 

At zero temperature, H is equal to the logarithm of the spe- 
cific enthalpy 



H = In 



e + P 



n B m B 



(8) 



where m B is the baryon mass. As a basic test of the numerical 
method of integration, we have checked, in the case of tabu- 
lated EOS at T = 0, the agreement between H calculated nu- 
merically and the analytical result (8). Due to the finite number 
of points in the tabular EOS the agreement is not complete, but 
the discrepancy was always found to be smaller than 1 % in all 
tested cases (even with strong differential rotation). 

3. Equation of state and thermodynamics 

Most of previous studies of structure of rotating neutron stars 
have only been performed with polytropes, or at most zero tem- 
perature EOS. This is not enough for the purpose of studying 
PNS evolution, since thermal and chemical effects are relevant 
and need a careful treatment. Only in a few works (Goussard 
et al. 1997, 1998; Strobel et al. 1999) a realistic, finite tem- 
perature EOS has been used. The works by Goussard et al. 
(1997, 1998) used a Skyrme-like potential model (Lattimer & 
Swesty 1991), while the more recent work by Strobel et al. 
(1999) used a Thomas-Fermi model with momentum and den- 
sity dependent effective nucleon-nucleon interaction developed 
by Myers & Swiatecky (1990). 

In this paper, we employ an EOS in the framework of rela- 
tivistic field theory, in which the interaction between nucleons 
is mediated by the exchange of three meson fields. The cou- 
pling constants are chosen to reproduce some bulk properties of 
nuclear matter such as the nuclear saturation density, symme- 
try energy or nuclear compressibility. The formalism is widely 
used, and to avoid unnecessary repetition we refer to the com- 
plete review by Prakash et al. (1997) and references therein. 
The particular choice of the coupling constants is detailed in 
Pons et al. (1999). 
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Figure 1 shows profiles from ID simulations of a 
M B = 1.6 M & spherical star from Pons et al. (1999). In pre- 
vious works (Goussard et al. 1997, 1998; Strobel et al. 1999) 
the evolution of the star was approximated by choosing either 
simplified profiles of constant entropy or temperature, or the 
combination of a constant low entropy core (s — 1 - 2) sur- 
rounded by a high entropy mantle (s = 4 - 5) to mimic the sit- 
uation in the early Kelvin-Helmholtz stage. It was shown that 
the choice of the critical density that separates the hot envelope 
from the inner core affected substantially the results. To reduce 
the uncertainties associated to the choice of this critical density 
and the fact that PNSs are not isothermal or isentropic config- 
urations, we have chosen to use the thermodynamical structure 
from ID simulations as a better approximation to reality, and 
assume that at each time the temperature and lepton fraction 
have the same dependence on the baryon density as the non- 
rotating models, T = T[n B ] and Y L = Y L [n B ]. More explicitly: 
given the results from ID simulations, the profiles of density, 
temperature, lepton fraction, entropy and pressure are tabulated 
at each time step, and these tables are taken as effective one- 
parameter EOSs where the independent variable is the baryon 
density (n B ) and the rest of variables are functions of only n B . 

Notice that, as we explained in the previous section, this 
(pseudo-barotropicity) assumption is a sufficient condition to 
ensure that the equation of motion has a first integral, making 
affordable the computation of the models. Of course, to know 
the exact profiles consistently, one should solve the neutrino 
transport equations in 2D along the whole evolution, but this 
requires a large technical and computational effort that is out 
of the scope of this paper. 

4. Choice of the initial models 

One of the first issues that must be addressed to study the early 
evolution of neutron stars is the initial model. In Pons et al. 
(1999) it was shown that for spherical models the overall evo- 
lution is not very sensitive to modifications of the initial ther- 
modynamical profiles, while the total mass of the star is the pa- 
rameter that mostly affects the subsequent evolution. It seems 
therefore natural to chose a canonical star with a baryon mass 
1 .6 Mq, that corresponds to a gravitational mass of the old, cold 
configuration of about 1.44 M . For completeness, we have 
also studied a model with a baryonic mass of 1.2 M . Notice 
that one cannot use the effective barotropes obtained for a PNSs 
with a given mass in the calculations of sequences with a differ- 
ent mass, because the thermodynamic properties are different. 

More difficult is to guess the particular rotation properties 
of a newly born neutron star. Newly born neutron stars are ex- 
pected to rotate non rigidly, as opposite to old neutron stars in 
which the various viscous mechanisms had enough time to act. 
Even in the case that the iron core of the Supernova progenitor 
was rigidly rotating (a quite reasonable assumption, since it is 
an old massive star), the process of collapse would generate a 
significant amount of differential rotation. Indeed, recent cal- 
culations of stellar evolution indicate that the iron core of mas- 
sive stars is rotating almost rigidly (Heger et al. 2000, 2003) 
and this is the approximation used in the most recent simu- 
lations of rotating stellar core collapse (Miiller et al. 2003). 



Based on this works, we have decided to study a few mod- 
els from Dimmelmeier et al. (2002, DFM hereafter) and to run 
some more simulations of similar models of stellar core col- 
lapse varying the total initial angular momentum. Our models 
correspond to those labeled by Al in DFM, i.e., the iron core of 
the progenitor is almost rigidly rotating. In Figure 2, we show 
profiles of the angular velocity (upper panel) and density (lower 
panel) of the inner 150 km, as a function of the equatorial ra- 
dius, immediately after the PNS is formed. The different lines 
correspond to models with a different amount of initial angu- 
lar momentum, with values of the ratio between rotational to 
gravitational potential energy (|r/W|) of 0.9% (models B3 in 
DFM), 0.5% (B2), 0.25% (Bl) and 0.05%, respectively. We in- 
cluded the last lower value because it is more realistic (Heger 
et al. 2003; Miiller et al. 2003) than those studied in DFM. In 
the figure we also show with dots the simple fit 



Q[r,#] = Q[r sin)?] 



+ r A sin 



(9) 



which corresponds to the Newtonian limit of O consistent with 
Eq. (4). In this limiting case, we now see that Ro is the value 
of the radius at which, in the equatorial plane (§~=7r/2), the an- 
gular velocity is half its limiting central value, Q c . Even if this 
remark is valid only in the Newtonian case, it is nicely veri- 
fied for relativistic stars as well. For all four models we have 
taken fl 2 = 50 km 2 , while Q e takes the values of 4500, 3500, 
2600, and 1300 rad/s, from top to bottom. The agreement be- 
tween the simple law (9) and the results from simulations is 
quite acceptable. In all the models the spatial scale of variation 
of the angular velocity (Ro) is the same because the density pro- 
files were very similar (see Figure 2, bottom). Varying the EOS 
might lead to a different density distribution, in such a way that 
one should expect smaller Ro for softer EOSs (more compact 
PNS), and larger Ro for stiffer EOSs. We must remark again 
that our assumption of having stationary axisymmetric solu- 
tions implied that the function F on the r.h.s. of the equation of 
stationary motion (2) is only a function of the angular velocity. 
In the Newtonian case, this is equivalent to say that the angular 
velocity depends only on the distance to the axis (rsini?). In 
the relativistic case, there are small corrections but the overall 
distribution of angular velocity is nearly cylindrical, as can be 
seen in Fig. 3. 

In reality, immediately after core collapse and bounce, there 
is no reason to expect that the distribution of angular velocity 
corresponds to that of stationary equilibrium, and it may take 
several rotation periods for the star to relax to some station- 
ary solution. A more detailed analysis of all these issues is out 
of the scope of this paper, and requires a parametric study of 
rotating core collapse simulations. However, for our purpose, 
which is just to understand in which range of parameters one 
should move when using Eq. (4), the conclusion is that a rea- 
sonable realistic model with differential rotation should be con- 
sistent with Rq of the order of 10 km. Considering that there 
might be a number of processes that could lead to rigid rotation 
on a dynamical timescale, such as shear or magneto-rotational 
instabilities, we will study both cases, with substantial differ- 
ential rotation and with rigid rotation. Notice that previous 
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works (Goussard et al. 1998) on differentially rotating PNSs 
have focused on the parameter space region corresponding to 
Rq « 1 km, which seems too low according to our results. 

A final remark concerns the geometry of the PNSs. 
Although for progenitors with moderate angular momentum, 
conservation of angular momentum leads to a characteristic 
profile of angular velocity that decreases with increasing dis- 
tance to the rotation axis, in the case of very fast rotation, or un- 
realistic large differential rotation, the usual oblate shape of the 
PNS can become toroidal, with a maximum in density off axis. 
In this case the profile of the angular velocity can exhibit also 
a maximum at some distance from the axis. We will discuss 
below for which models a toroidal shape is found, but unfortu- 
nately our numerical algorithm does not allow us to compute 
equilibrium configurations with zero density at the origin. 

5. Results 

We have considered a number of thermodynamical profiles cor- 
responding to snapshots of the early evolution of PNSs, and we 
have calculated the rotational properties in two different cases: 
rigidly rotating and differentially rotating PNSs. The left part of 
Table 1 shows some properties (central baryonic density, grav- 
itational mass, and radius) of the non rotating models with a 
baryonic mass of 1.6 M Q . The models are labeled according 
to the evolutionary time (in s after birth) in the non-rotating 
case. Notice that for rotating PNSs the evolution time label- 
ing the models is only an indication, since in a self-consistent 
simulation one would expect faster evolution. The reason of 
this faster evolution is that the lower densities, and consistently 
lower temperatures because the entropy is conserved, result in 
larger neutrino mean free paths. Therefore the neutrino trans- 
port would be more effective. 

In order to investigate the influence of the low density 
(< 10 11 g/cm 3 ) EOS, we have also substituted the zero tem- 
perature BPS equation of state (Baym et al. 1971) by the LS 
(Lattimer & Swesty 1991) EOS in some models at early times. 
The results for the models with the LS EOS are shown between 
parenthesis in the first four rows of Table 1 . At later times, af- 
ter the mantle loses most of its thermal and lepton content and 
contracts, the influence of the low density EOS becomes less 
important and both EOSs give essentially the same result, be- 
cause the amount of mass and the thickness of the low density 
layer is rather small. The left hand side of Table 2 shows the 
equivalent information but for a non rotating NS with a total 
baryonic mass of 1.2 M Q . 

5.1. Rigidly rotating PNSs 

In the right hand side of Tables 1 and 2 we summarize the prop- 
erties of maximally rotating PNSs following the same thermo- 
dynamical sequence as in the non-rotating case. Focusing at 
early times (f < 0.5 s) we see that the central densities of non- 
rotating and maximally rotating PNSs are not terribly different. 
The reason is that these objects are only slightly bounded and 
they cannot rotate very fast, reaching the mass shedding limit at 
angular velocities of hundreds of Hz. Therefore, the rotational 
energy is not large and the models are not very different from 



the non-rotating case. This can also be seen in the values of 
the quantity \T/W\, which is about 0.03-0.05 for the BPS EOS 
and 0.06-0.07 when using the LS EOS at low density. The dif- 
ferences between these two EOSs come from the fact that the 
BPS is at zero temperature while the LS is at finite temperature. 
Due to the more extended mantle, the total angular momentum 
is systematically large, of about 30%, for the models with finite 
temperature at p < 10 11 g/cm 3 . 

After about 0.5-1.0 seconds, the mantle has contracted and 
lost its thermal and lepton content, so that the PNS becomes 
more compact and bound. This can be deduced by looking at 
the central density, which exhibits a sudden increase at the end 
of this first epoch. For the non-rotating models it happens be- 
tween 0.5 and 1.0 seconds, while for the maximally rotating 
models it seems to be delayed until some time between 1.0 
and 2.0 seconds, due to the role of the centrifugal force. For 
these more compact configurations, the limiting angular veloc- 
ities are much larger (3000-6000 rad/s), and their central densi- 
ties are considerably lower than those of the non-rotating PNSs. 
The ratio between polar and equatorial radii of maximally ro- 
tating PNSs turns out to be rather insensitive to the thermo- 
dynamics, being about 0.65 for early models and decreasing 
to 0.55 for late configurations. With respect to the \T/W\ ra- 
tio, its maximum value is about 0. 1 at the end of the evolution, 
still below the threshold to undergo dynamical instabilities. As 
a rule of thumb, we found that the maximum angular veloc- 
ity can be estimated with good accuracy with the simple fit 
Q K ~ 0.58 ^GM/R\ valid for both masses (1.6 and 1.2 M ) 
and all thermodynamical profiles studied. Despite our models 
are not close to the maximum mass for a given EOS, this result 
is in agreement with the empirical formula derived by Haensel 
& Zdunik (1989) relating the maximum rotation frequency to 
the mass and radius of the maximum mass static configuration. 
It seems that this empirical relation is valid in general. 

The evolution of the angular momentum deserves a sepa- 
rate discussion. The first result is that including thermal effects 
at low density (LS) allows for significantly larger total angu- 
lar momentum. However, this discrepancy disappears rapidly, 
and the maximum angular momentum decreases quite fast and 
monotonously with time, which has evolutionary implications. 
Whether or not neutrino transport can redistribute and take 
away a substantial fraction of the initial angular momentum is 
unclear and needs of multidimensional transport simulations to 
be clarified. In general, it is assumed that the angular momen- 
tum is approximately constant during the Kelvin-Helmholtz 
epoch, and it only varies on timescales larger than the diffu- 
sion timescale. If this was the case, PNSs born with maximal 
angular momentum cannot radiate it via neutrino diffusion with 
the same efficiency as the binding energy, and this may result 
in significant mass loss after reaching the mass shedding limit. 
Other possibility is that PNSs are not born maximally rotating, 
the angular momentum is (almost) conserved in a timescale of 
seconds, and the star speeds up as it contracts and evolves ther- 
modynamically. We have simulated this scenario by taking a 
sequence of configurations with constant angular momentum 
( J - 1 .5 GMq/c), always with fixed baryonic mass. Results are 
shown in Figure 4. It can be seen how a model that initially had 
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a \T/W\ ratio of only 0.02 speeds up in a timescale of a sec- 
ond and then slowly approaches the mass shedding limit in a 
few seconds. We expect that fast (but not maximally) rotating 
PNSs might reach the threshold of some types of instabilities 
(CFS, bar-mode, ...) on a diffusion timescale, say within 10 s 
after birth. It will depend on the efficiency of neutrinos to carry 
away angular momentum, which is currently under investiga- 
tion. 

The above discussion was focused on results for the canon- 
ical NS with a baryonic mass of 1.6 M Q . For completeness we 
also give results for a model with a lower baryonic mass of 1 .2 
M in Table 2. The results are qualitatively similar to those of 
Table 1, with only some quantitative differences. Accordingly 
to the lower mass, the maximum angular velocities and \T/W\ 
ratio are also lower. 

5.2. Differentially rotating PNSs 

As discussed in the previous section we will restrict to realis- 
tic values of Rq in the range from 10 km (significant differen- 
tial rotation) to 50 km (almost rigid rotation). A different sce- 
nario would be a hot NS resulting from the merger of two neu- 
tron stars. In this case, the distribution of angular momentum 
is more complex and depends not only on the initial total an- 
gular momentum, but also on a number of different parameters 
as the mass ratio or the relative orientation of the spins of the 
two stars. Lower values of Ro might be reached in this latter 
scenario, but above all, the rotation profile would probably be 
very different from Eq. (4), which is in agreement with a sta- 
tionary configuration, the case of a merging being far from a 
stationary situation. 

As before, the baryonic mass of the star is conserved and 
fixed to 1.6 M G for all models. Within the range of parameters 
mentioned above, we have calculated a large number of mod- 
els of differentially rotating PNSs with either a fixed R of 10, 
20, 50, and 100 km or varying Rq according to the circumfer- 
ential equatorial radius of the star (Rq = 0.5 R eq ). In Figures 
5 and 6 we show iso-enthalpy surfaces for two representative 
models at an evolutionary time of t = 1 s. Despite both config- 
urations have similar equatorial radii 30 km), they are qual- 
itatively very different. In Figure 5 the parameter Rq — 100 km, 
i.e., the model has a very low degree of differential rotation. 
Even very close to the mass shedding limit (0.995 of the max- 
imum angular velocity) the star shows a regular oblate profile. 
On the contrary, the model in Figure 6 corresponds to the case 
of significant differential rotation (Rq = 10 km). With strong 
differential rotation we have found that above some critical an- 
gular velocity (950 rad/s, for this model) the shape of the star 
becomes toroidal, with a maximum density off axis. 

The most relevant results are presented in a compact form 
in Figures 7 and 8, where we show the total angular momentum 
and the rotation parameter \T/W\ as a function of the central 
angular velocity £2 C , for different values of Ro- For compari- 
son, we have also included the results for rigid rotation (solid 
curve), together with those for Ro = 50 km (dash-dotted line), 
Ro = 20 km (thick dashes), Ro = R eq /2 (thin dashes) and 



Rq = 10 km (dots). For Ro > 100 km, the results are so close to 
the rigid rotation case that we do not show them in the figures. 

The thermodynamical profiles for each panel are those cor- 
responding to evolutionary times of the non-rotating models of 
t = 0.5 (top), 1.0 (middle) and 10 s (bottom panel) after for- 
mation. From the results of the figure we can see that for Rq 
of the order or larger than 50 km the effects of differential ro- 
tation are unimportant. Only for Rq < 20 km and Rq = 0.5R eq 
the differences are significant. As discussed in 4, results from 
axisymmetric core collapse seem to indicate that the scale in 
which one should expect variations of the angular velocity is of 
about 10 km, therefore the last two cases may be looked at as 
the most realistic, with the caveat of being constrained to the 
particular rotation law we used. 

In the model corresponding to t — 1 s it is visible a sudden 
increase of the angular momentum at a critical value of the an- 
gular velocity, different for each value of Rq. This is an effect of 
the phase transition from bulk matter to nuclei, and it is asso- 
ciated to the point when the angular velocity becomes so large 
that the central density is below nuclear saturation density (po). 
This transition is not observed at earlier times because matter 
at low density is still very hot, and the thermal effects suppress 
or smooth out the kink. At later times, the star is more com- 
pact and the central density is always above po, even close to 
the mass shedding limit. For the models with strong differen- 
tial rotation, we had to stop our calculations before reaching 
the mass shedding limit because of numerical limitations asso- 
ciated to the special toroidal geometry that arises (see Fig. 6). 
These cases are marked by an asterisk in Fig. 7 and 8. In prin- 
ciple, there could be configurations with slightly larger angular 
momentum. 

An expected result is that the central angular velocity can 
be up to a factor 5 to 10 larger in the case of differential ro- 
tation, but the maximum angular momentum is only about a 
fifty percent larger than for rigid rotation. More important is the 
variation in the maximum of \T/W\. While for rigidly rotating 
stars it is about 0.1, for PNSs with a large degree of differen- 
tial rotation it can be as high as 0.2. The fact that differential 
rotation allows for the existence of equilibrium models at large 
\T/W\ have also been found by Baumgarte et al. (2000) and 
Yoshida et al. (2002) for polytropes in the relativistic Cowling 
approximation using the same rotation law. A maximum value 
of \T/W\ = 0.2 is not enough to reach the dynamical instability 
threshold (« 0.27), but it is sufficient to allow either the secu- 
lar, gravitational wave driven, instability (Chandrasekhar 1970; 
Friedman & Schutz 1978) that happens at 0.1 - 0.15 or the re- 
cently proposed low \T/W\ dynamical instabilities (Shibata et 
al. 2002, 2003; Watts et al. 2003). Notice that despite our cal- 
culations cannot be carried up to the mass shedding limit for 
some models, we are not far from it and the numbers are close 
to the maximum value. This can be seen in Figure 9, when we 
show the polar to equatorial axis ratio as a function of Q e . For 
rigid rotation, or moderate differential rotation the calculations 
end at the Keplerian limit, while for strong differential rotation 
the calculation ends because of numerical problems (marked 
with asterisks). The slope of the curves is becoming more and 
more vertical, and the axis ratio is as low as 0.4, smaller than 
the minimum in the rigid case (0.6). This is an indication that 
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Table 1. Properties of non-rotating PNSs and rigidly rotating PNSs at the limiting frequency, for a fixed baryonic mass M B = 1.6 Mq. The 
models are labeled by the evolutionary time of non-rotating PNSs at which the thermodynamical profile was calculated. The entries in the 
table are: central baryon number density, n c ; gravitational mass, M a ; circumferential equatorial radius, R eq ; Kepler frequency, f2 K ; angular 
momentum, J; axis ratio r (polar to equatorial) and the rotation parameter \T/W\ (rotational energy on absolute value of gravitational energy). 
For 4 time steps at early times we give results using the LS EOS at low densities for comparison (first four rows, between parenthesis). 
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Table 2. Properties of non-rotating PNSs and rigidly rotating PNSs at the limiting frequency, for a fixed baryonic mass M B = 1.2 Mq. The 
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the shape of the star is changing fast and only slightly larger 
angular velocities can be reached before the star starts loosing 
mass. 

The next question one can pose is whether or not dynami- 
cal instabilities arise in a real case. As for rigidly rotating stars, 
we now impose that the angular momentum is fixed along the 
evolution, as well as the baryonic mass, and study how the dif- 
ferent variables evolve under such constraints. Figure 10 shows 
the temporal evolution of the central angular velocity (top) and 
of the \T/W\ ratio (bottom) of a PNS differentially rotating with 
Ro = 10 km and a fixed angular momentum of J = 1.5 GM^/c. 
As in the case of rigid rotation, \T/W\ increases as the star 
contracts and looses its binding energy in a neutrino diffusion 
timescale of about 10 s. Thus, if there are no significant losses 
of angular momentum, the various types of dynamical insta- 
bilities associated to critical values of the rotation parameter 
may arise not at the very beginning, but several seconds after 
formation. 

As a consistency check of our assumptions and many sim- 
plifications, one can make a simple estimate of the neutrino 
luminosity considering that the variation in the gravitational 
mass of the PNS is considered to be equal to the luminosity, i.e., 
L v = -^jf-- The result is shown in Figure 11. Both the or- 
der of magnitude and the exponential decay of the luminosity 
in a timescale of 10 seconds are consistent with the results of 
detailed simulations with neutrino transport in spherical sym- 
metry (Pons et al. 1999). It is interesting that the difference 
between the two models studied with Ro = 10 km (solid) and 
Ro = 20 km (dashes) is only visible after 10 seconds of evolu- 
tion, when the stars become more compact. Therefore, impos- 
ing conservation of angular momentum and baryonic mass, and 
building a sequence keeping fixed Ro is not violating any basic 
physical law, such as conservation of energy, and can be taken 
as a qualitative approach to the real case. 

6. Conclusions 

We have approached the problem of the evolution of rotating 
protoneutron stars by constructing evolutionary sequences of 
axisymmetric stationary configurations in General Relativity. 
The thermodynamical structure and evolution have been taken 
and extrapolated from simulations in spherical symmetry that 
included neutrino transport. Although this is a crude simplifi- 
cation, it already gives an interesting insight about how the dif- 
ferent relevant quantities can evolve as the rotating PNS loses 
its lepton content and its excess binding energy, and contracts. 
Moreover, we have found that the luminosity estimates are not 
terribly different from what one expects. 

A special effort has been made to understand in which 
space of parameters we should be in a realistic case. The 
biggest uncertainty concerns the rotation law that PNS have at 
birth. By analyzing results from simplified core collapse sim- 
ulations, it seems that a typical scale for variations of the an- 
gular velocity is about 10 km, and that conservation of angu- 
lar momentum during the collapse of a stellar core (initially 
rigidly rotating) does not seem to allow for angular velocities 
varying in a length-scale shorter than a few km. Less is known 
about the angular distribution, except that the most recent simu- 



lations show the presence of important meridional currents and 
some turbulent motion. For simplicity, we restricted ourselves 
to the stationary case. Stationarity implies a quasi-cylindrical 
distribution (with deviations due to relativistic corrections) of 
the angular velocity. This stage only can be reached after sev- 
eral dynamical and rotation periods, after the PNS had time 
to relax. Therefore one must be aware that the first 0.5 s are 
probably far from stationarity, but after that evolution proceeds 
in a quasi-stationary way, except for low velocity convective 
motions. From our study of quasi-stationary sequences we can 
draw a few interesting qualitative results. 

i) For rigidly rotating stars, the mass shedding limit is well 
approximated by the simple law Cl K « 0.58 ^GM/R 3 , despite 
complications in the EOSs. This empirical formula is valid for 
stars of different masses within a 5% accuracy for the EOSs we 
considered. We believe that this result is general and would be 
valid for other EOSs or stars with different masses. 

ii) For each instant in the evolution, stars with strong dif- 
ferential rotation can have 5 to 10 times larger central angular 
velocities, and accommodate about a fifty percent more angular 
momentum. The maximum specific angular momentum J/M 
varies between (1 - 2) GM B /c * (0.5 - 1) x 10 16 cm 2 /s, de- 
pending on the degree of differential rotation. 

iii) When R Q < 20 km there is a critical value of the angular 
velocity above which the shape of the star becomes toroidal, 
with a maximum density off axis. We found that this config- 
urations are stationary solutions of the hydrodynamics equa- 
tions, but it is unclear under which circumstances this geometry 
is stable against perturbations. This issue deserves a separate 
study. 

iv) The maximum value of \T/W\ obtained in the case of 
differential rotation is about 0.2, while for rigid rotation this 
is a 0.11. Thus differentially rotating PNSs might undergo the 
CFS instability, which arises at « 0.14 and, in any case, the 
recently discussed low \T/W\ instability (Shibata et al. 2002, 
2003; Watts et al. 2003) is plausible to happen. 

v) More interestingly, we found several situations in which, 
even if the initial model is not close enough to the critical value 
of \T/W\, as the star contracts in a neutrino diffusion timescale 
of 5-10 s, it speeds up entering the window of instability. An 
observational evidence of this effect could be a temporal shift 
of a few seconds between the neutrino luminosity peak and a 
gravitational wave burst in the event of a galactic Supernova. 
Ultimately, this depends on the initial amount of angular mo- 
mentum, which is approximately equal to the angular momen- 
tum of the iron core of the progenitor. Recent stellar evolu- 
tion calculations suggest that the specific angular momentum 
of the inner 1.7 M Q of a 15 M Q star can be as high as 3 x 10 16 
cm 2 /s if magnetic braking is neglected, or 10 15 cm 2 /s if mag- 
netic torques are included in the evolution (Heger et al. 2003). 
This corresponds to J/M in the range J/M = (0.2 - 6) GM G /c. 
If the angular momentum happens to be in the upper region 
( J/M > 2), centrifugal forces would stop the collapse before 
the PNS is formed. Intermediate values (J/M =1) may result 
in the formation of a rapidly rotating PNS that enters the insta- 
bility region several seconds after birth. If magnetic braking is 
very effective during the evolution of a massive star, J/M < 0.5 
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and the PNS will be formed after collapse without reaching ex- 
treme values of the angular velocities and \T/W\. 

The next natural step to improve this work is to include the 
possible mechanisms to transport angular momentum between 
the different layers of the star, that may involve neutrino trans- 
port, turbulent transport, magnetic fields, neutrino viscosity, 
convective motion and/or angular momentum losses by grav- 
itational wave emission. Unless the star is born with almost 
maximal angular velocities, since the different instabilities can 
arise in a timescale of seconds some of all of these dissipative 
mechanisms can modify our current vision of PNS evolution. 
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Fig. 2. Equatorial profiles of angular velocity (top), and density (bot- 
tom) of PNSs from axisymmetric simulations of stellar core collapse 
(DFM). The four models correspond to a different amount of angu- 
lar momentum of the iron core, namely, \T/W\ = 0.9% (solid), 0.5% 
(dashes), 0.25% (dash-dot), and 0.05% (dash-3 dots). The dotted lines 
on the upper panel are fits to the simple law (9) with Rl = 50 km 2 . 
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Fig. 1. Temporal variations of the temperature (top), entropy (middle) 
and lepton fraction (bottom) as functions of the baryon density, for 
evolutionary times corresponding to t = 0.1,0.5, 1 and 8 seconds after 
formation of the PNS. Results are for a M B = 1.6M G star in spherical 
symmetry. 
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Fig. 3. Surfaces of constant angular velocity (projected in a transversal 
plane) of a differentially rotating star with R = 100 km, for t = 1 s, 
and rotating with CI = 354 rad/s. (Q./Q K = 0.995). The surface is 
marked by the thick solid line. Because of relativistic effects, the sta- 
tionary solution does not show exact cylindrical symmetry. 
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Fig. 4. Temporal evolution of the angular velocity £l/Q K and the ro- 
tation parameter for a sequence of rigidly rotating PNSs with 
constant angular momentum (1.5GMq/c) and a fixed baryonic mass 
of M B = 1.6 Mq. 
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Fig. 5. Iso-enthalpy lines of a differentially rotating star with 
R = 100 km, for t = 1 s, and rotating with Q. = 354 rad/s. (C1/Q K = 
0.995). For low differential rotation the structure of the star still looks 
normal. 
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Fig. 6. Iso-enthalpy lines of a differentially rotating star with 
Ro = 10 km, for t = 1 s, and rotating with f2 = 1570 rad/s. The 
toroidal shape is clearly visible for all models differentially rotating 
with R smaller than a certain threshold. 



Fig. 7. Total angular momentum as a function of the central angular 
velocity for PNSs with fixed baryonic mass of 1.6 M Q at evolutionary 
times of 0.5 (top), 1.0 (middle) and 10 s (bottom), and for different 
values of R : R = oo (solid), 50 km (dash-dot), 20 km (thick dashed), 
R eq /2 (thin dashed), 10 km (dots). 
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Fig. 9. Axis ratio (polar to equatorial) as a function of the central an- 
gular velocity for a PNS with a baryonic mass of 1.6 M Q at t = 10 s, 
and for different values of R : R Q = oo (solid), 50 km (dash-dot), 20 
km (thick dashed), R eq /2 (thin dashed), 10 km (dots). 
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Fig. 8. Rotation parameter, \T/W\, as a function of the central angular 
velocity for PNS with fixed baryonic mass of 1.6 M Q at evolutionary 
times of 0.5 (top), 1.0 (middle) and 10 s (bottom), and for different 
values of _/? : = 00 (solid), 50 km (dash-dot), 20 km (thick dashed), 
R eq /2 (thin dashed), 10 km (dots). 
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Fig. 11. Estimated neutrino luminosity (-dMo/dt) as a function 
of time, consistent with the evolution at fixed baryonic mass of 
M B = 1.6 Mq and angular momentum of 1.5 GM^/c. The two models 
are differentially rotating with Ro = 10 km (solid line) and Ro = 20 km 
(dashed line). The order of magnitude and the exponential decay are 
similar to the luminosities obtained in simulations with neutrino trans- 
port for non-rotating stars. 



Fig. 10. Temporal evolution of the central angular velocity Q c and of 
the rotation parameter \T/W\ for two sequences of differentially rotat- 
ing PNSs with Rq = 10 km (solid line) and Ro = 20 km (dashed line). 
In both cases we fix the total angular momentum to 1.5 GM^/c and 
the baryonic mass to M B = 1.6 Mq. 



